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Abstract. We model the transport of cosmic ray nuclei in the Galaxy by means 
of a new numerical code. Differently from previous numerical models we account 
for a generic spatial distribution of the diffusion coefficient. We found that in the 
case of radially uniform diffusion, the main secondary/primary ratios (B/C, N/0 and 
sub-Fe/Fe) and the modulated antiproton spectrum match consistently the available 
observations. Convection and re-acceleration do not seem to be required in the energy 
range we consider: 1 < -E < 10^ GeV/nucleon. We generalize these results accounting 
for radial dependence of the diffusion coefficient, which is assumed to trace that of the 
cosmic ray sources. While this does not affect the prediction of secondary/primary 
ratios, the simulated longitude profile of the diffuse 7-ray emission is significantly 
different from the uniform case and may agree with EGRET measurements without 
invoking ad hoc assumptions on the galactic gas density distribution. 



1. Introduction 



Understanding the origin and propagation of Cosmic Rays (CR) in the Galaxy is an 
intriguing puzzle which requires the combination of many different observations over 
a wide range of energies. While simplified propagation models, most noticeably the 
leaky-box and the modified slab models, succeeded interpreting the main features of 
the CR nuclear composition and energy spectra for 1 < ^ < 100 GeV/ifl, more 
realistic diffusion models are called for to provide a comprehensive description of 
multi-channel observations (including heavy nuclei, electrons, 7-rays and antimatter 
particles) accounting for the growing amount of available astrophysical data (see [Ij for 
a comprehensive review). 



J This unit of measure corresponds to the energy E competing on average to each nucleon within a 
nucleus with A nucleons and energy A x E. 
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Two main approaches have been developed so far: analytical (or semi-analytical) 
diffusion models (see e.g. [2] and ref.s therein), which solve the CR transport equation 
by assuming simplified distributions for the sources and the interstellar gas, and fully 
numerical diffusion models. Well known recent realizations of those two approaches are 
respectively the two-zone model [31 H] and the GALPROP code [HI El [7] . In the case 
of GALPROP, the adoption of realistic gas and radiation field distributions allows to 
model also the spectrum and angular distribution of the 7-ray secondary emission. 

Although these models represent a significant step forward with respect to previous 
analyses, they still perform a number of simplifications with respect to a more realistic 
physical scenario. Most noticeably, they assume diffusion to be statistically isotropic 
and homogenous, i.e. they adopt the same diffusion coefficient value (or at most two 
values in the two-zone model) all over the propagation volume. 

However, such assumptions may not always be justified, as diffusion coefficients 
generally depend on the regular magnetic field orientation and on the ratio between 
the regular and chaotic magnetic field energy densities. Although these quantities are 
poorly known, several observations and theoretical arguments suggest that they are far 
from being spatially homogeneous in the Galaxy (see sec. [2] for more details). This may 
have relevant consequences for the CR spatial distribution in the Galaxy, for the angular 
distribution of the secondary 7-ray and neutrino emissions [8] and to interpret the CR 
anisotropy. 

In order to be able to test some of those effects, as well as to verify previous 
results which have been derived in the literature under more conventional conditions, 
we developed a new numerical code, DRAGON (Diffusion of cosmic RAys in 
Galaxy modelizatiON). DRAGON is especially designed to account for a spatially in- 
homogeneous and an-isotropic diffusion coefficient. In its present version it allow to 
model CR nuclei transport at energies -Emin ^ 1 GeV/n as well as the secondary 7-ray 
and neutrino emission produced by their interaction with the interstellar medium (ISM). 
In this work we disregard CR convection and re-acceleration (i.e. we work in a plain 
diffusion (PD) regime) and show that most relevant measurements can be reproduced 
under these conditions. Above £^mm we expect that no other physical input than source 
spectra, diffusion and fragmentation processes can determine secondary/primary ratios, 
hence a comparison of our prediction with experimental data should allow to fix the slope 
delta of the diffusion coefficient for some assumed slope of the CR injection spectra (see 
[9] for a detailed discussion about this issue). 

In our analysis we will mainly refer to measurements of the secondary/primary flux 
ratios of several nuclear species (the most relevant are B/C, N/0 and sub-Fe/Fe) and 
the antiproton and p/p spectra, performed by several satellite and balloon experiments. 

In order to test our code, we firstly study the conventional case of a uniform diffusion 
coefficient. Afterwards, we will analyze the previously unconsidered case in which D 
grows exponentially with the distance from the Galactic Plane (CP) and traces the 
radial distribution of supernova remnants (SNR). 

This paper is structured as follows: in sec. [2] we give several theoretical and 
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observational motivations for considering the effects of in-homogeneous diffusion on 
the CR distribution. Moreover, we explain why it is worth studying perpendicular, 
rather than isotropic diffusion. In sec. [3] we describe our numerical model and our main 
assumptions on several variables entering the CR transport equations. In particular, 
we describe our assumed CR source and target gas distributions, and the nuclear cross 
section models we exploit. Then, we describe our main results on the physics of CRs, by 
studying secondary/primary ratios and individual species spectra in sec. HI We will show 
that the CR spatial distribution can be significantly affected by in-homogeneous diffusion 
and discuss how this can have relevant consequence for the solution of the CR gradient 
problem, as we already argued in [8]. Although to simulate a detailed map of the 7-ray 
diffuse emission, suitable to comparison with the one the FERMI observatory will soon 
provide jlO], is beyond the aims of this work, in sec. [5] we will model the longitudinal 
profile of this emission along the CP showing that under reasonable condition it matches 
EGRET measurements. Finally, in sec. |6] we draw our conclusions. 



2. Theoretical and observational motivations for in-homogeneous diffusion 
models 

Charged particles diffuse in chaotic magnetic fields due to their scattering onto hydro- 
magnetic fluctuations. The presence of a regular component of the magnetic field, which 
is the case in the Milky Way, is expected to break isotropy so that spatial diffusion has 
to be described in terms of a diffusion tensor L)jj(x). According to [11] this can be 
conveniently decomposed as 

Ai(x) = (Di(x) - D||(x)) BiBj + D\\{^)5ij + DA{^)tijkBk , (1) 

where Bi are the components of the regular magnetic field versor. The symmetric 
components D\\ and are the diffusion coefficients along and perpendicularly to the 
regular field Bq, while Da is the antisymmetric (Hall) diffusion coefficient which accounts 
for the drift due to the interplay of Bq and CR density gradient. Since Da is relevant 
only at very high energies (i? > 1 PeV, see e.g. [H [12]) we will disregard it in the 
following. 

Since diffusion is related to magnetic processes, diffusion coefficients depend on the 
particle rigidity p = p{E)/Ze. Moreover, in general D\\ and D± depend differently on p 
and on the strength of hydro- magnetic fluctuations. In the quasi-linear theory (QLT) 

D\\{^,p)c^^vrL{p)V~\p) (2) 



J^. 5B{kf 

where ri(p) = p/Bo is the Larmor radius and V{p) = — „ is the integral of the 

normalized power spectrum of the turbulent hydromagnetic fluctuation over modes with 
wavenumber k > fcmin(p) = 27rr^^(p). A power-law behavior B{kY oc k''^ is generally 
assumed, with 7 = 5/3 (3/2) for Kolmogorov (Kraichnan) turbulence spectrum implying 
D|| oc p^/^ (1/2) In QLT the perpendicular diffusion coefficient is 

D^(x,p)~D||(x,p)P(p)2«l, (3) 
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meaning that diffusion takes place mainly along the regular magnetic field lines. 

Although QLT may not be applicable to the conditions presents in the ISM, 
more realistic computations |2j confirmed that expectation founding D± ~ 0.1 D\\. 
MonteCarlo simulations of particle propagation in turbulent fields [131 [I2l [H] also found 
a similar result (although computation time limits allowed to test it only at energies 
above 100 TeV). 

What is most relevant here, however, is the different behavior of D\\ and D± as a 
function of the turbulent power. Simulations of propagation in strongly turbulent fields 
agree with QLT predicting D\\ {D±) decreasing (increasing) when V{p) increases. It 
should be noted that if, as it is generally assumed, the CR source distribution can be 
approximated to be cylindrically symmetric, and the regular field to be purely azimuthal 
B = (0, 0, Bff,), parallel diffusion plays no rol^. 

Clearly, under this approximation and in the absence of an a priori criterion to fix 
the normalization and energy dependence of the diffusion coefficients, the substitution 
of an isotropic diffusion coefficient with D± would produce no physical effects. This 
conclusion is no more true, however, if the homogeneous diffusion approximation is 
relaxed and one tries to correlate spatial variations of the relevant diffusion coefficients to 
those of the hydro-magnetic fluctuation energy density, as Dy and D± have an opposite 
behavior as functions of 'P(p). 

Observationally, very little is known about the spatial distribution of hydro- 
magnetic fluctuations in the Galaxy. There are, however, evidences both for a longitude 
[T6] and latitude dependence of the fluctuation power [171 HE] . Also from a theoretical 
point of view, 'P(p) is quite unlikely to be uniform as fluctuations are expected to be 
correlated, via particle-wave resonant scattering, to CRs which, in turn, are correlated 
to the non-uniform source distribution. 

A radial variation of the diffusion coefficient may have relevant consequences on 
the CR spatial distribution in the galactic disk. In [8] some of us already pointed out 
that in-homogeneous diffusion may help reconciling the discrepancy between the rather 
smooth diffuse 7-ray longitude profile observed by EGRET [I9] with the quite steep 
SNR (the most likely CR sources) radial distribution {CR gradient problem). That can 
be understood as a back-reaction effect: a larger CR density nearby sources induces a 
larger V{p), hence a larger D±, which in turn implies a faster CR diffusion out of those 
regions (note that the effect would be opposite for Dy). In sec. |5]we will discuss this 
effect in more details and show its possible relevance for the 7-ray angular distribution. 

Concerning the vertical profile of the diffusion coefficient, we assume here 

D{p,r,z) = D{p,r) exp{\z\/zt} . (4) 

This behaviour, which was also adopted in [21j, naturally leads to a more regular profile 
of the CR density at large \z\ (see figure [11). We note, however, that as far as stable 

§ This conclusion is not expected to change significantly if a possible spiral shape and a tiny dipole 
component of the regular magnetic fields are accounted for (though more complex scenarios have been 
considered [T5]V 
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Figure 1. The proton flux vertical profile at 1 GeV obtained with DRAGON assuming 
a uniform diffusion coefficient (black, dashed line) is compared with that obtained 
adopting the exponential profile in equation |4] for zt — 4: kpc (blue, continuos line). In 
both cases D is normalized so to reproduce the B/C (see section |4|). 



secondary nuclei are concerned, replacing an exponentially vertically growing D± with 
a uniform one has almost no effects, as expected because spallation takes place mainly 
in the thin Galactic disk where the CR density is only marginally affected by the choice 
between these two options. 



3. Description of the model 

It is well known since the pioneering work of Ginzburg and Syrovatskii [22], that in 
absence of continuos energy losses, re-acceleration and convection, the diffusive transport 
of stable nuclei in the ISM is described by the following equation 

^ + V.iD.VN,)= (5) 

Qi{Ek, r,z) -cp ragas(r, z) ai^{Ek)Ni + ^c(3 ng^{r, z) ajiNj . 

j>i 

where Ek = {E—mAj/A {E is the total energy of a nucleus with mass 771 A — ^ X 'HT.jmcleon ) 
is the kinetic energy per nucleon E^, which is constant during propagation and is 
practically conserved in fragmentation reactions. (3 is the velocity of the nucleus in 
units of the speed of light c. CTj is the total inelastic cross section onto the ISM gas 
with density r;,gas(r, 2) and aij is the production cross-section of a nuclear species j 
by the fragmentation of the i-th one. We start the spallation routine from A = 64. 
We disregard ionization and Coulomb energy losses as they are known to be negligible 
above ~ 1 GeV/n. Indeed, we verified with GALPROP that their contribution to the 
predicted secondary /primary and antiproton spectra is only a few percent above this 
energy. 
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We solve equation numerically in the stationary limit dNi/dt = using a 
Cranck-Nicholson scheme [231 [211 [8] with open boundary conditions. This corresponds 
to free escape of CRs at the outer limit of the Galaxy, defined by -Rmax = 20 kpc and 
Zraa.x- While -Rmax IS fixed, Zmax IS Set to z^aay: = 2 X 2:^ (cf. cq. Hj) to avoid bordcr effects. 

We describe below our assumptions for the terms appearing in eq. ([5]). 

Spatial diffusion The l.h.s. of eq. ((51) describes spatial diffusion. As we discussed 
above, we assume cylindrical symmetry and that the regular magnetic field is 
azimuthally oriented (Bq = B^{r,z) 4>)- Under these conditions CR diffusion out 
of the Galaxy takes place only perpendicularly to Bq. Therefore in the following it 
is understood that D represents in fact the perpendicular diffusion coefficient D±. 
The dependence of D on the particle rigidity p is (see e.g. [25] ) 



The function /(r) describes a possible radial dependence of D. We define it to be 
unity at Sun position (r = r©) so that Dq corresponds to the local value of the 
diffusion coefficient at the reference rigidity po = 3 GV. 

CR sources For the source term we assume the general form 



imposing Qi{Ek,rQ, Zq) = 1. While the exact form of fs{r,z) has little influence 
on the charged secondary species spectra, it is very important in shaping the 7-ray 
angular distribution. We assume fs{r, z) to trace the SNR distribution as modeled 
in [26] on the basis of pulsar and progenitor star surveys [8]. In the galactic disk 
such a distribution is similar to that adopted in [27], but shows an excess in the 
Galactic Bulge due to the contribution of type-la SNe, not accounted for in [27]. 
Both distributions are significantly more peaked than those empirically determined 
[23, E] by matching the 7-ray longitude profile measured by EGRET [TU] . 
The injection abundances are tuned so that the propagated spectra of primary 
and secondary (or their ratio) species fit the observed ones (see below). Even 
though our code allows to consider different power-law indexes for the different 
nuclear species, in this work we only consider the same ai = a for all species, 
when not differently stated. For each value of 5 in eq. (El) the source spectral slope 
a is fixed by the requirement that at high energy ^ 100 GeV/n, at which 
spallation processes are almost irrelevant, the equality a + 6 = 2.7 is satisfieclii], 
in order to match the observed slope. It is understood that DRAGON (similarly 
to GALPROP) cannot account for the local conditions (e.g. enhanced turbulence 
strength and gas density) in the vicinity of the sources, because its limited spatial 
resolution does not allow to investigate lengths < 40 pc (while the typical source 

II In this regime, the theoretical expectation for the observed flux $ on Earth is ^{E) w Q[E) / D{E) ~ 




(6) 




(7) 
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size is c/ < 10 pc). These effects, however, are unhkely to produce a significant 
effect on the CR spectrum, as the mean fraction of the CR path-length spent in 
those regions is quite smaU. 

Nuclear cross sections As in GALPROP (see |7j and references therein) the 
spallation cross sections and the spallation network are based on a compilation 
of experimental data and semi-empirical energy dependent interpolation formulas 
as provided e.g. in [291 EOl l3T] . 

Target gas The IS gas is composed mainly by molecular, atomic and ionized hydrogen 
(respectively, H2, HI and HII). Although more realistic distributions are known, for 
r > 2 kpc we adopt the same distributions as in GALPROP, for essentially two 
reasons. First of all, since CRs propagate for million years in the Galaxy, in the 
stationary limit they just probe a smoothed, mean gas distribution. Secondly, we 
can have a more direct comparison with GALPROP results. 

However, in the central region of the Galaxy, where GALPROP assumes an 
interpolated density, we use the the H2 and HI distributions as modeled in [32] . 
While the flux and composition of charged CR reaching the Earth are not sensitive 
to the central gas distribution, this choice allows us to better model the 7-ray 
emission in the Galactic Centre (GC) region, as we will discuss in more details in 
sec. m Following [31] we take the He/H numerical fraction in the ISM to be 0.11. 
We neglect heavier nuclear species. 

4. Testing DRAGON: the case of a radially uniform diffusion coefficient 

In order to test our code, we ran it under similar conditions to those already considered 
in the literature. In this section we show the results we obtained assuming that the 
diffusion coefficient does not depend on the galactocentric radius r. As we mentioned in 
sec. m the adoption of an exponential vertical profile for D does not affect significantly 
the results presented in this section with respect to the case of isotropic and uniform 
diffusion mostly considered in literature. Indeed, passing from a spatially uniform D to 
the profile described by eq. ([6]) only amounts to a small re-scaling of Dq. 

In the following, every label indicating a nucleus refers in fact to the sum of all its 
isotopes, unless otherwise stated. 

4.1. The B/C ratio 

The Boron to Carbon ratio (B/C) is one of the most useful tracers of CR propagation in 
the Galaxy. In fact, since Boron is entirely secondary, its observed abundance strongly 
depends on the residence time of primary CRs in the Galaxy. Moreover, measurements 
of Boron and Carbon fluxes are better than those of other secondary /primary ratios, 
and the B production cross sections from its main primaries (^^C and ^^O) are known 
better than for other secondary nuclides. 
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4- 1.1. Fixing free parameters Once the spatial distributions of the CR sources and 
the ISM gas have been chosen, the main parameters determining the B/C in a PD 
model are the C/0 and N/0 injection ratios and the quantities S, Dq and Zt in eq. (E]). 
As it was already shown in several papers, secondary/primary ratios for stable species 
depend on the ratio Dq/zi (which will be always expressed in units of 10^^ cm^ s~^ kpc~^ 
throughout this paper) rather than on the two parameters separately. 

While primary /primary ratios are usually disregarded in the literature, as they do 
not give direct relevant information on CR propagation, we use them to fix the C/0 and 
N/Oljll injection ratios, while we fix the abundances of primaries heavier than oxygen by 
requiring that they match the observed abundances in CRs at ~ 1 — 10 GeV/n. 

To this aim, we define two different . We compute the former (which we label 
X{c/o,N/o}) t>y comparing our predictions for the C/0 and N/0 modulated ratios to 
experimental data over the energy range of our interest. The latter (which we label 
^{Do/zt &}) computed comparing our predicted modulated B/C ratio to the observed 
one. Solar modulation is taken into account here in the "force-field" approximation [33] 
using a modulation potential of magnitude $ = 500 MV. In order to study potential 
energy dependent effects we consider two different minimum kinetic energies per nucleon 
E^™ for comparison to data: 1 GeV/n and 2 GeV/n. The low statistical significance of 
the data set above this energy prevents us from going further up the energy scale. 

For each pair of values {Do/zt, S), we determine the X{c/o,n/o} distribution in the 
space (C/0, N/0) scanning over a wide range of C/0 and N/0 injection ratios. For 
the set of parameters that minimizes X{c/on/o}' compute xfoo/zt s} ^^"^ repeat 
this procedure for several values of {DQ/zt,S). Finally, we analyze the distribution of 
X{Do/zt s} obtain our best fit values for {Dq/zi, 6) with the appropriate confidence 
regions. 

Thus, this strategy allows us to fix best values of the C/0 and N/0 injection ratios 
and to consistently determine the best propagation parameters that will be used as our 
best model for the analysis of antiproton and 7-ray fluxes. 

We notice here that this procedure, which corresponds essentially to split the whole 
4:—D parameter space into two separate ones, is physically motivated by the weak 
dependence of primary /primary ratios on {Do/zt,S). 

4-1.2. Experimental data So far the best B/C measurements above 1 GeV/n have 
been provided by the HEAO-3 [35] and CRN [36] experiments in the range 1 < Ek < 
30 GeV/n and 70 GeV/n < Ek < 1.1 TeV/n. New data should be released soon by 
the CREAM [3^ and TRACER [2H] experiments significantly improving the available 
statistics at high energy. Here we use only HEAO-3, CRN data. For consistency, we 
take also C/0 and N/0 data from the same experiments. 

^ Note that N — ^^N + ^^N is a combination of primary and secondary nuclides 
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4- 1.3. Results We show in figures 2(a) and 2(b) the main results of our procedure. 
In fig. 2(a) the distribution of X{Do/zt,5} -E™" = 1 GeV/n is shown, together 



with confidence regions at 68%, 95% and 99% Confidence Level (CL). Our best-fit 
values for {Do/ Zt, 5,C/0,N /O) if Ef"" = 1 GeV/n are (0.55, 0.57, 0.79, 0.044). The 
projection of this point in the {Do/zt,S) plane is highlighted by the cross in figure 



2(a) Remarkably, the best-fit value for S favors a Kraichnan turbulence spectrum, 
rather than a Kolmogorov one. Changing the minimum energy -E™" from 1 GeV/n to 
2 GeV/n indeed produces no relevant effect. In particular, the best-fit values for D^/zt 
and 6 are not moved (see fig. 2(b)). It is interesting to notice that the particular value 
of 6 = 0.57 we obtain is consistent with findings of other authors (see [1] and references 
therein). The best-fit C/0 and N/0 injection ratios (0.79,0.044) should be compared 
with the solar system ones [34j 0.76 and 0.11 respectively. 

We do not include the sub-Fe/Fe (sub-Fe = V -|- Ti -|- Sc) ratio in our statistical 
analysis because of the large uncertainties on the knowledge of the spallation cross 
sections for heavy elements. However, we found that we consistently obtain a reasonable 
match of experimental data also for this observable. In order to improve the fit to this 
ratio a careful fine tuning of nuclear cross section parametrizations seems to be needed. 
In fig. [3] we show the B/C, C/0 and N/0 ratios as obtained with our best-fit model, 
and assuming 2:4 = 4 kpc. 

A comment is in order here: the particular observables we are considering are 
not sensitive to Dq and Zt independently. A possible way to estimate Zt is offered by 
unstable/stable ratios (e.g. ^°Be/^Be), which are known to probe the vertical height of 
the Galaxy p] . Unfortunately, the best experimental data for this particular ratio have 
been obtained at energies < 100 MeV/n, while only 2 experimental points with large 
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Figure 3. In these panels we show our best fit for the B/C, C/0 and N/0 compared 
with HEAO-3 [Sg (red diamonds), CRN [35] (green, triangle) and ATIC-II [Mi (blue) 
experimental data (though the latter are not used in our statistical analysis). Continuos 
curves: local interstellar (LIS); dashed lines: top of atmosphere (TOA) ($ = 500 MV). 



errors are available at 1 GeV/n [IQ]. Since our code does not allow us to have reliable 
predictions down to few hundreds MeV/n, because we do not account for energy losses, 
it is impossible for us to draw any firm conclusion about our best value for Zt. However, 
by considering only the ^°Be/^Be points around 1 GeV/n we infer that Zf should lie 
between 3 and 5 kpc, in agreement with previous findings [S]. 

4-2. Antiprotons 

Most antiprotons reaching the Earth are expected to be a product of CR hadronic 
collisions with the IS gas. Their measured spectra provide, therefore, valuable 
information on CR propagation which are complementary to that coming from secondary 
nuclei (see e.g. plHHlHi]). 

The main processes responsible for p production are p Pgas, p Hcgas, He Pgas and 
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Figure 4. The p/p ratio (left) and antiproton absolute spectrum (right) are compared 
with BESS 95+97 i48j, BESS98 gS], CAPRICE [50], and AMS [ST] experimental data. 
The shadowed regions correspond to models matching the B/C data within la. LIS: 
light yellow band, between continuos lines; modulated (<I> = 550 MV): grey band, 
between dashed lines. 



He Hegas, plus a negligible contribution from other nuclei. Similarly to [131111] we use the 
p production cross-section calculated using the parametrization given in Tan & Ng [45j . 
We account for the contribution of heavier nuclei in the CRs and the ISM by using the 
effective correction function determined by Simon et al. [16] with the MonteCarlo model 
DTUNUC. Inelastic scattering, annihilation and tertiary p (antiprotons which have been 
inelastically scattered) are treated as in [S]. For the local interstellar spectrum (LIS) of 
primary nuclei we adopt $p = 1.6 x 10^ {Ek/1 GeV)"^'^'^ (m^ s sr GeV)~^ as measured 
by BESS during the 1998 flight [17] by accounting for a solar modulation potential 
$ = 550 MV in the "force- free" approximation [33]. 

We use DRAGON to simulate the primary proton distribution in the Galaxy and the 
LIS of secondary antiprotons. Normalization is imposed by requiring that the simulated 
proton LIS coincides with $p(rQ, 0)obJil- In figures 4(b) and 4(a) we compare our results 



with the experimental data released by BESS for the periods 1995-97 [IE] and 1998 [19] 
in the energy interval 1-4 GeV, and by CAPRICE (1998) [50] in the range 3-49 GeV. 
All these data refer to a period of low solar activity (the minimum was in 1997) and 
same positive phase of the solar cycle, with a mean value of the modulation potential 
for the period 1995-1998 of $ = 550 MV p]. Hence we will also use $ = 550 MV to 
obtain p modulated spectra. 

The shaded regions correspond to the uncertainty on the antiproton flux due to 
the uncertainty on the propagation parameters, and are constructed using values of 
{Do/zt, S) within the la region of fig. 2(a) The figures shown in this section are derived 
using Zt = 4 kpc, but we checked that, similarly to the B/C case, different choices of this 
parameter do not affect the antiproton spectrum provided that Do/zt is kept constant. 

+ Only the absolute p flux is dependent on such normalization. 
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In fig. 4(a) we also show the LIS and TOA energy behavior of the p/p ratio obtained 



with the parameters corresponding to minimum of the distribution shown in fig. 2 (a' 



It is evident that the models which fit the B/C data within la are also compatible with 
the antiproton measurements. A statistically poor excess of the predictions of our best- 
fit model respect to the BESS data, which was also found in [H], is probably not 
significant due to the large systematic uncertainties. At higher energies, we have a 
small tension between our predictions and the highest energy CAPRICE data. A better 
agreement may be found if preliminary PAMELA [52] results [53] will be confirmed. 



5. Radial dependent diffusion and the 7-ray longitude distribution 

In this section we model the secondary 7-ray emission originated, via tt" decay, by the 
interaction of the hadronic component of CRs with the IS gas. Along the GP, where the 
gas column density is higher, this process is expected to give the dominant contribution 
to the total diffuse emission above the GeV. At the energies of our interest a simple 
scaling model for the differential production cross section can still reliably be used. In 
this regime the energy spectrum of secondary 7's is a power law with the same slope as 
the primary nuclei (only protons and He nuclei give a significant contribution). 

The main gas (target) components are the molecular (H2) and atomic (HI) 
hydrogen, and He atoms. The contribution of ionized hydrogen is almost irrelevant 
in the GP. For r > 2 kpc, we adopt the same HI spatial distribution as [5J. For the H2 
we assume 

nH2(r, z) = eo(r) Xco(r) exp |- ln2(z - 2;o)V2:ft(r)} , (8) 

where e(r) is the CO (a widely used H2 tracer) volume emissivity, ZQ{r) and Zh{r) are 
the midplane displacement and scale heights respectively, and Xco('") is the CO - H2 
conversion factor. All these quantities, with the exception of Xco('"), are the same as in 
[5^l5|[6] for r > 2 kpc, while for smaller radii we adopt the Ferriere et al. model [32]. The 
adoption of Ferriere's model for the molecular and atomic hydrogen for r < 2 kpc allows 
us to avoid the interpolation of the 7-ray fiux profile in the GC region and to reproduce 
naturally the peaked emission observed by EGRET toward the GC as we already pointed 
out in [55]. For the 11% He fraction we adopt the same spatial distribution as for the 
HI. 



5.1. The CR gradient problem 

To simulate a detailed map of the 7-ray diffuse emission of the Galaxy (even only for the 
hadronic component) is beyond the aims of this work. Rather, the main issue we want to 
address here is the so called CR gradient problem. This originates from the well known 
discrepancy between the theoretical fiux profile obtained by assuming SNRs to be the 
sources of galactic CRs and that inferred from EGRET 7-ray diffuse observations [19] . 
Under mild assumptions on the distribution of the galactic gas, it was found [28], [5] that 
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the inferred CR radial profile should be much flatter than the theoretically expected 
one. 



5.1.1. A former proposed solution and a test for DRAGON A possible way out was 
suggested in [27] in terms of a radially variable Xqq. While in i6j this quantity was 
assumed to be uniform (Xco = 1-8 x 10^° cm~^/(K kms"^))El, in [27] it was taken to 
increase gradually by more than one order of magnitude from 4 x 10^^ at r = 2 kpc to 
1 X 10^^ at r > 10 kpc. However, while the growth of this parameter with r is suggested 
both by theoretical arguments and observations of external galaxies, its actual behavior 
is rather uncertain so that in [27J it had to be tuned into 5 steps to match EGRET 
observations. 

To test our code against possible failures in reproducing the 7-ray longitude profile, 
we try to reproduce the results of |27j. We adopt the same XQo{r) which was used in 
[27] and a CR model giving the best-fit of the B/C in the case of a radially uniform 
diffusion coefficient. We use 5 = 0.57 (see sec. 14.1.31) but our results do not change 
appreciably by using any value in the interval 0.45 — 0.65. In fig. 6(a) we compare our 
results with EGRET measurements along the GP for A < < 10 GeV [56] . 

We reasonably reproduce both the normalization and the main features of the 
observed longitude profile. Smaller structures may only be reproduced using a detailed 
3—1^ model of gas distribution which we are planning to do in a forthcoming paper. For 
comparison, in the same figure we also show the emission profile which we would obtain 
using a constant Xco('") = 1.8 x lO^'' for r > 2 kpc. 



5.2. An alternative solution of the GR gradient problem 

As an alternative possibility we explore the case in which the diffusion coefficient traces 
the radial dependence of the SNR distribution as we motivated in sec. [21 According 
to the arguments explained in the same section we expect the CR radial profile to be 
smoothed with respect to the one obtained in the case of constant diffusion coefficient. 
Hence, we expect to be able to fit EGRET longitude profile without fine tuning the 
parameter Xco- Indeed, this is what we find. 

We assume a constant Xqq = 1.8 x 10^° for r > 2 kpc, while in the bulge 
(r < 2 kpc), where physical conditions are much different from the outer disk, we 
take Xco = 0.5 x 10^° [32] . For the diffusion coefficient, we assume that the function 
/(r), as defined in eq. ([6]), is 

f{r) = fs{r,Or . (9) 

The function fs{r, 0) describing the radial distribution of the Galactic SNRs is taken by 
[26] and is the same as in [8J. 

A self-consistent model should be developed in order to rigourously determine 
/(r) and the primary CR density Xj(r) accounting for the mutual influence of CR 

* For clarity, in the following we will drop units in quoting values of Xco ■ They are always understood 
to be cm~^/(K knis^^). 
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and hydromagnetic fluctuations. Although such an analysis is beyond the aims of this 
work, a simple estimate can be derived under reasonable physical assumptions. In 
particular, we assume that, as observed locally and expected by theoretical reasons, 
energy equipartition holds between CR and magnetic field fluctuations. Under this 
hypothesis, and restricting our analysis to the galactic disc, we have that the fluctuation 
strength V is proportional to the CR density A^(r, 0), and, since equations ([2]) and ([3]) 
imply Dx{r) oc V, we obtain -D_L(r) oc A^(r, 0) (we further simplify the problem by 
assuming that D_l is constant in z, since, as we also pointed out in the previous section, 
our observables do not depend strongly on the vertical behaviour of Dx). Assuming 
now that the source density is of the form fs{r,z) ~ g{z)F{r), it is straightforward to 
show that, in the limit in which radial diffusion can be neglected with respect to the 
vertical one (or, in other words, that the diffusion halo radius is much larger than is 
vertical thickness), the stationary limit of the spatial diffusion equation takes the form 

■ — — = Kg{z)F{r), where i^' is a constant. This equation can be easily 
solved by assuming N{r, z) = N{r, 0)H{z), from which we obtain 

D^oc N{r,0) = F{rf-\ (10) 

implying r = 0.5. 

Of course, this rough estimate cannot replace a more detailed self-consistent 
numerical computation. It is intriguing, however, that, besides its pedagogical utility, 
this simple derivation predicts a CR radial distribution which is quite close to what 
is needed to solve the gradient problem. In figure [5] we show the radial profile of the 
proton differential flux at 1 GeV for three different values of the parameter r: 0,75, 0.5 
and (the latter corresponds to a radially uniform D). It is clear that for r > 0.5 a 
significant flattening of the CR spatial distribution with respect to that of sources has 
to be expected. 

In order to verify if such distributions are compatible with CR measurments, we 
run DRAGON under the same conditions discussed in sec. H] but adopting a radially 
dependent D as specified in ([9]). We find that, as long as r < 1, we are still able to 
obtain a good fit of secondary/primary ratios and antiproton data. In particular, we 
find that for r ~ 0.5 this happens with the same propagation parameters which allow to 
match observations in the case r = 0, in spite of the fact that those values correspond 
to quite different CR radial distributions. 

The degeneracy of the different CR density profiles with respect to CR 
measurements could be removed by means of 7-ray observations if combined with 



independent astrophysical measurements of the Xco radial dependence. In fig. 6(b) 



we show the simulated 7-ray longitude profile as obtained using r = 0.5 and Xco taken 
to be uniform for r > 2 kpc. Clearly this is quite smoother than the profile obtained 
under the same conditions but a radially uniform D [r = 0). A slightly better fit of data 
can be obtained using r = 0.75, probably indicating that also radial diffusion should 
be taken into account respect to our semplified analytical model . The comparison of 
figures 6(a) and 6(b) displays the degeneracy between the radial dependence of Xco and 
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Figure 5. Proton differential flux at E = 1 GeV for three different choices of the 
parameter t setting the radial dependence of the diffusion coefficient on the SNR 
distribution (see eq. (O). r = (radially uniform D), blue, continuous curve; r = 0.5, 
red, dash-dotted; r = 0.75, violet, dotted. In all cases Zt = 4: kpc and the D 
normalization giving the best fit to B/C data is chosen. 



that of the diffusion coefficient. It is intriguing, however, that the EGRET profile may 
be explained in terms of in-homogeneous diffusion with no need to invoke a tuning of the 
Xqo parameter. While the actual 7-ray longitude profile is likely to be determined by a 
combination of the radial dependence of both those quantities, it should be clear from 
our results that the effect pointed out in this work should be taken into account when 
interpreting observations of 7-ray diffusion emission of the Milky Way and of external 
galaxies. 

6. Conclusions 

We modeled CR transport in the Galaxy assuming a plain diffusion model. We study 
two main cases: in the first one the diffusion coefficient D is assumed to be uniform 
along the Galactic Plane, while in the second case we consider, for the fist time, a D{r) 
which traces the radial profile of SNRs (which we assume to be the CR sources). 

The fact that we can consistently reproduce the observed antiproton spectra and the 
main secondary /primary nuclear ratios for E >1 GeV/n makes us quite confident of the 
validity of our approach. By using only the B/C, C/0 and N/0 data we found that the 
preferred range (la) of values of the slope of diffusion coefficient is 0.43-0.65. The best 
fit value is 5 ^ 0.57. This is in agreement with findings of other authors. A Kolmogorov 
spectrum is disfavored and re-acceleration seems to be unnecessary to interpret data 
above 1 — 2 GeV/n. Forthcoming experiments like CREAM |37j and TRACER [38] for 
what concerns nuclei and PAMELA ||52] and AMS [5l] for antiprotons may soon allow 
to strengthen this conclusion by improving both statistics and quality of data. 

While in both cases (uniform and radially dependent D) we obtain substantially the 
same successful predictions for what concerns nuclei and antiprotons reaching the Earth, 
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Figure 6. Our predictions for the longitudinal profiles of the 7-ray hadronic emission 
integrated for |6| < 1° are compared with EGRET measurements [19]. Left panel: 
radially uniform D with Xco — 1-8 x 10^° cm^^/(K kms~^) for r > 2 kpc (long 
dashed curve), and Xqo as in [27 for r > 2 kpc (continuous curve). Right panel: 
D{r) tracing the SNR distribution with t = 0.5. In both cases Zt = A kpc and D 
normalization is chosen to best-fit the B/C data. 



the corresponding CR primary spatial distributions in the Galaxy can be considerably 
different. This may have a number of interesting effects, including a possible role in 
the solution of the problem which plain diffusion models face predicting a too high CR 
anisotropy above 100 TeV. 

We focused here on the effects on the expected secondary 7-ray diffuse emission. 
We showed that the longitude distribution of that emission can be significantly affected 
by in-homogeneous diffusion. In [8] we already noticed that the effect goes in the right 
direction to provide a viable solution of the CR gradient problem. Here we provide 
further theoretical arguments in favour of a radially dependent diffusion coefficient and 
succeed reproducing EGRET observations for A < E < 10 GeV and |6| < 1° for a 
choice of the relevant diffusion parameters which allow to match the B/C and the 
antiproton spectrum. The extension of our predictions to larger latitudes would require 
to implement in DRAGON electron propagation (and losses) and more detailed gas and 
radiation distributions which we plan to do in a forthcoming paper. 
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